Returns on Investment to Herbicides: Example of Spatial Economic Surplus Approach

Author

Maxwell Mkondiwa

1 Introduction

The conventional approach to valuing agricultural research is the economic surplus approach (Alston et al 1995). This approach is commonly applied to country wide research programs rather than disaggregated to district level. In this note, we demonstrate how to use causal machine learning heterogeneous treatment effects to make return on investment analyses that are sufficiently disaggregated.

2 Steps in Return on Investment Economic Surplus Model

2.1 Step 1: Baseline area, production, and yield

Code
library(rio)

APY_data <- import("Data_1997_23_Cleaned.xlsx", sheet = "All_data")

APY_data$Production_Tonnes <- as.numeric(APY_data$Production_Tonnes)

APY_data_wheat_2015_19 <- subset(APY_data, APY_data$Season == "Rabi" & APY_data$Crop == "Wheat" & APY_data$Year %in% c("2015-16", "2016-17", "2017-18", "2018-19", "2019-20"))

APY_data_wheat_2015_19$Season <- NULL
APY_data_wheat_2015_19$Crop <- NULL
APY_data_wheat_2015_19$Season <- NULL
library(data.table)
APY_data_wheat_2015_19 <- data.table(APY_data_wheat_2015_19)
APY_data_wheat_2015_19_ave <- APY_data_wheat_2015_19[, lapply(.SD, base::mean, na.rm = TRUE), by = .(Districts)]

APY_data_wheat_2015_19_ave <- APY_data_wheat_2015_19[, .(Area_Hectare = base::mean(Area_Hectare), Production_Tonnes = base::mean(Production_Tonnes), Yield_Tonnes_Hectare = base::mean(Yield_Tonnes_Hectare)), by = .(Districts)]

# Download all the data
library(reactable)
library(htmltools)
library(fontawesome)

htmltools::browsable(
    tagList(
        tags$button(
            tagList(fontawesome::fa("download"), "Download as CSV"),
            onclick = "Reactable.downloadDataCSV('APY_data_wheat_2015_19_ave', 'APY_data_wheat_2015_19_ave.csv')"
        ),
        reactable(
            APY_data_wheat_2015_19_ave,
            searchable = TRUE,
            defaultPageSize = 38,
            elementId = "APY_data_wheat_2015_19_ave"
        )
    )
)
Code
# library(ggplot2)
# ggplot() +
#     geom_line(data = subset(APY_data, APY_data$Crop == "Wheat"), aes(x = Year, y = Area_Hectare, group = Districts, linetype = Districts, color = Districts), linewidth = 1.5) +
#     labs(x = "Year", y = "Area") +
#     theme_bw() 


# Having the baseline we need to create this from 2020 to 2030

library(tidyverse)
APY_data_wheat_2019_30 <- APY_data_wheat_2015_19_ave %>%
    mutate(year = list(2019:2030)) %>%
    separate_rows(year, sep = ",", convert = T)

# Delete Kishanganji due to lack of data
APY_data_wheat_2019_30 <- subset(APY_data_wheat_2019_30, !(APY_data_wheat_2019_30$Districts == "KISHANGANJ"))

2.2 Step 2: Adoption data and projections

Code
# District level Cumulative Adoption -------------------------------------------
library(rio)
DistrictCumulativeAdoption=import("DistrictCumulativeAdoption.csv") 

# Manual
library(nlme)
DistrictCumulativeAdoption100=nlsList(cumulative_percent~100/(1 + exp(-(a + b*herbi_1st_use))),start=list(a = 1, b = 0),
                                      data=groupedData(cumulative_percent~herbi_1st_use|A.q103_district,DistrictCumulativeAdoption),na.action=na.omit)
summary(DistrictCumulativeAdoption100)
Call:
  Model: cumulative_percent ~ 100/(1 + exp(-(a + b * herbi_1st_use))) | A.q103_district 
   Data: groupedData(cumulative_percent ~ herbi_1st_use | A.q103_district,      DistrictCumulativeAdoption) 

Coefficients:
   a 
                Estimate Std. Error    t value     Pr(>|t|)
Madhubani      -368.5966   82.71001  -4.456493 1.027323e-03
Sheohar        -440.8519   69.53585  -6.339922 1.196379e-05
Samastipur     -497.3282   85.30907  -5.829722 4.128779e-04
Saharsa        -822.5614  108.00732  -7.615793 1.719296e-05
Gopalganj      -737.3012   99.81949  -7.386345 9.473423e-04
Gaya           -465.0329   52.77353  -8.811860 4.595092e-07
Banka          -611.1577   67.42975  -9.063621 3.680051e-07
EastChamparan  -555.9987   64.48144  -8.622616 1.265154e-06
Saran          -801.4401  107.21386  -7.475153 1.303373e-03
Arah           -550.8303   71.37689  -7.717208 3.420563e-04
Kaimur         -615.2736   85.27350  -7.215297 2.174164e-04
Sitamarhi      -605.5380   82.36076  -7.352263 1.293976e-04
Bhagalpur      -549.7067   54.12638 -10.155985 6.352263e-08
Siwan          -527.0180   73.48044  -7.172221 9.610372e-05
Patna          -623.5362   82.43138  -7.564307 3.007824e-04
Madhepura      -868.6877   99.10018  -8.765753 1.487608e-06
Muzaffarpur    -798.2624   92.32364  -8.646348 1.547847e-05
Munger         -766.7486   95.12688  -8.060273 7.326273e-06
Araria         -316.5343   51.04528  -6.201050 1.498404e-02
Supaul        -1132.2963  143.64343  -7.882688 5.506730e-05
Khagaria       -741.4620   82.52171  -8.985054 2.718708e-06
Begusarai      -676.2813   68.29952  -9.901698 1.427091e-07
Katihar        -840.8264  109.69662  -7.665016 5.133371e-04
Rohtas         -552.4232   53.60328 -10.305772 2.646961e-09
Nalanda        -552.0948   51.53724 -10.712542 2.859514e-08
Lakhisarai     -660.6143   65.34005 -10.110405 9.762755e-07
Jehanabad      -572.6040   56.34564 -10.162349 2.685139e-06
Buxar          -870.8090   98.04538  -8.881693 6.783745e-06
Nawada         -626.5890   73.07958  -8.574064 3.771172e-04
Arwal          -609.8696   69.55790  -8.767798 7.598773e-05
Vaishali       -561.0634   67.36747  -8.328403 7.363271e-03
Aurangabad     -340.2565   31.86115 -10.679356 2.510768e-05
Darbhanga     -3624.7742  818.85530  -4.426636 1.803145e-02
Kishanganj            NA         NA         NA           NA
WestChamparan  -919.0901  129.75365  -7.083347 3.671990e-02
   b 
               Estimate Std. Error   t value     Pr(>|t|)
Madhubani     0.1824767 0.04103458  4.446901 1.035827e-03
Sheohar       0.2185873 0.03450839  6.334322 1.203434e-05
Samastipur    0.2465753 0.04232892  5.825221 4.146130e-04
Saharsa       0.4077597 0.05355931  7.613238 1.723119e-05
Gopalganj     0.3657399 0.04952904  7.384353 9.486854e-04
Gaya          0.2308192 0.02621096  8.806209 4.620421e-07
Banka         0.3031446 0.03345784  9.060494 3.691977e-07
EastChamparan 0.2759868 0.03201262  8.621187 1.266763e-06
Saran         0.3974994 0.05318910  7.473325 1.305022e-03
Arah          0.2733708 0.03542745  7.716357 3.423079e-04
Kaimur        0.3053386 0.04231911  7.215146 2.174447e-04
Sitamarhi     0.3006063 0.04089236  7.351162 1.294910e-04
Bhagalpur     0.2727610 0.02685720 10.155976 6.352316e-08
Siwan         0.2614366 0.03646276  7.169961 9.627491e-05
Patna         0.3095687 0.04092676  7.563969 3.008563e-04
Madhepura     0.4307406 0.04915224  8.763397 1.490672e-06
Muzaffarpur   0.3960573 0.04580312  8.646951 1.546978e-05
Munger        0.3805968 0.04722806  8.058700 7.338021e-06
Araria        0.1572961 0.02534405  6.206433 1.494838e-02
Supaul        0.5617378 0.07126618  7.882250 5.508221e-05
Khagaria      0.3678722 0.04094198  8.985207 2.718288e-06
Begusarai     0.3355991 0.03389792  9.900285 1.429328e-07
Katihar       0.4172890 0.05443688  7.665557 5.131419e-04
Rohtas        0.2743447 0.02662110 10.305535 2.647610e-09
Nalanda       0.2740957 0.02558715 10.712241 2.860414e-08
Lakhisarai    0.3279545 0.03243563 10.110932 9.758054e-07
Jehanabad     0.2841554 0.02795783 10.163713 2.681869e-06
Buxar         0.4323910 0.04867673  8.882908 6.777524e-06
Nawada        0.3110665 0.03626996  8.576424 3.765523e-04
Arwal         0.3028315 0.03452546  8.771251 7.577619e-05
Vaishali      0.2787775 0.03345523  8.332852 7.347176e-03
Aurangabad    0.1693484 0.01585108 10.683708 2.502638e-05
Darbhanga     1.7986934 0.40631770  4.426815 1.803072e-02
Kishanganj           NA         NA        NA           NA
WestChamparan 0.4563168 0.06439980  7.085687 3.668355e-02

Residual standard error: 5.446441 on 250 degrees of freedom
Code
alphaDistrictCumulativeAdoption100 <- coef(DistrictCumulativeAdoption100) 
write.csv(alphaDistrictCumulativeAdoption100, file = "alphaDistrictCumulativeAdoption100.csv")

# SSLOGIS
SSlogDistrictCumulativeAdoption100=nlsList(cumulative_percent~SSlogis(herbi_1st_use, 100, xmid, scale),start=list( xmid=1992, scale=5),
                                           data=groupedData(cumulative_percent~herbi_1st_use|A.q103_district,DistrictCumulativeAdoption),na.action=na.omit,pool=FALSE)

summary(SSlogDistrictCumulativeAdoption100)
Call:
  Model: cumulative_percent ~ SSlogis(herbi_1st_use, 100, xmid, scale) | A.q103_district 
   Data: groupedData(cumulative_percent ~ herbi_1st_use | A.q103_district,      DistrictCumulativeAdoption) 

Coefficients:
   xmid 
              Estimate Std. Error    t value     Pr(>|t|)
Madhubani     2019.965 0.62600442   3226.758 5.534587e-14
Sheohar       2016.823 0.29064217   6939.197 3.408871e-25
Samastipur    2016.942 0.38198596   5280.147 3.114780e-21
Saharsa       2017.270 0.18971836  10632.972 1.718670e-26
Gopalganj     2015.917 0.33389462   6037.584 1.393556e-21
Gaya          2014.706 0.28127301   7162.814 1.025929e-31
Banka         2016.060 0.22170394   9093.480 6.365014e-36
EastChamparan 2014.584 0.20750839   9708.446 1.419123e-29
Saran         2016.204 0.35377492   5699.116 1.969965e-21
Arah          2014.957 0.42340960   4758.883 4.257725e-27
Kaimur        2015.054 0.31172345   6464.235 5.599572e-25
Sitamarhi     2014.389 0.23554683   8551.966 4.149291e-19
Bhagalpur     2015.342 0.21881761   9210.144 5.603198e-36
Siwan         2015.854 0.28622508   7042.898 5.530895e-22
Patna         2014.209 0.31556520   6382.863 9.981837e-22
Madhepura     2016.731 0.15984861  12616.503 1.744640e-30
Muzaffarpur   2015.522 0.24927123   8085.660 3.446920e-32
Munger        2014.596 0.22275145   9044.142 1.257663e-32
Araria        2012.346 0.98121490   2050.872 2.556551e-10
Supaul        2015.703 0.11899314  16939.653 1.360759e-20
Khagaria      2015.543 0.24064125   8375.715 1.448376e-35
Begusarai     2015.146 0.22705890   8874.994 2.820863e-42
Katihar       2014.974 0.26861731   7501.281 3.788709e-22
Rohtas        2013.610 0.18035030  11164.996 3.737681e-40
Nalanda       2014.241 0.24769486   8131.945 8.055008e-42
Lakhisarai    2014.348 0.24608242   8185.662 1.822057e-35
Jehanabad     2015.109 0.33295856   6052.131 3.732649e-34
Buxar         2013.939 0.16696591  12061.976 7.109611e-27
Nawada        2014.325 0.43567904   4623.415 6.910769e-21
Arwal         2013.891 0.37380192   5387.589 1.577836e-27
Vaishali      2012.585 0.75791725   2655.415 1.437609e-16
Aurangabad    2009.211 0.68474963   2934.227 3.158098e-28
Darbhanga     2015.226 0.01337908 150625.217 4.226515e-06
Kishanganj          NA         NA         NA           NA
WestChamparan 2014.149 0.68307842   2948.635 7.937186e-14
   scale 
               Estimate Std. Error   t value     Pr(>|t|)
Madhubani     5.4801451 0.64234695  8.531441 1.035844e-03
Sheohar       4.5748162 0.41942225 10.907424 1.203451e-05
Samastipur    4.0555254 0.57710609  7.027348 4.146150e-04
Saharsa       2.4524195 0.23733817 10.333018 1.723180e-05
Gopalganj     2.7341816 0.45426114  6.018964 9.486895e-04
Gaya          4.3323885 0.34008558 12.739112 4.620479e-07
Banka         3.2987473 0.28183517 11.704527 3.692014e-07
EastChamparan 3.6233572 0.28186384 12.854991 1.266778e-06
Saran         2.5157257 0.44434601  5.661637 1.305020e-03
Arah          3.6580234 0.61488345  5.949133 3.423105e-04
Kaimur        3.2750447 0.46993566  6.969134 2.174477e-04
Sitamarhi     3.3265990 0.31398220 10.594865 1.294918e-04
Bhagalpur     3.6662154 0.26012456 14.094076 6.352441e-08
Siwan         3.8250014 0.41832460  9.143620 9.627511e-05
Patna         3.2302883 0.43344524  7.452587 3.008599e-04
Madhepura     2.3215746 0.18447540 12.584738 1.490692e-06
Muzaffarpur   2.5248860 0.30183702  8.365064 1.547023e-05
Munger        2.6274484 0.28659268  9.167884 7.338174e-06
Araria        6.3574427 1.25799511  5.053631 1.494840e-02
Supaul        1.7801852 0.14084622 12.639211 5.508187e-05
Khagaria      2.7183180 0.28830505  9.428618 2.718303e-06
Begusarai     2.9797434 0.27387486 10.879945 1.429352e-07
Katihar       2.3964023 0.35472053  6.755747 5.131420e-04
Rohtas        3.6450463 0.21169845 17.218106 2.647668e-09
Nalanda       3.6483484 0.29007825 12.577118 2.860487e-08
Lakhisarai    3.0491957 0.28917699 10.544393 9.758114e-07
Jehanabad     3.5191837 0.37269062  9.442641 2.681894e-06
Buxar         2.3127176 0.19456947 11.886333 6.777520e-06
Nawada        3.2147323 0.44943329  7.152858 3.765589e-04
Arwal         3.3021531 0.44588977  7.405761 7.577658e-05
Vaishali      3.5870959 0.82429109  4.351734 7.347191e-03
Aurangabad    5.9049921 0.74957546  7.877782 2.502655e-05
Darbhanga     0.5559592 0.01575021 35.298526 1.803048e-02
Kishanganj           NA         NA        NA           NA
WestChamparan 2.1914333 0.70991098  3.086913 3.668419e-02
Code
nlssummary=summary(SSlogDistrictCumulativeAdoption100)

xmid=nlssummary$coefficients[,,1]
colnames(xmid) <- paste("xmid", colnames(xmid), sep = "_")
scale=nlssummary$coefficients[,,2]
colnames(scale) <- paste("xmid", colnames(scale), sep = "_")

SSlogDistrictCumulativeAdoption100Results <-data.frame(coefficients(SSlogDistrictCumulativeAdoption100),xmid,scale)

library(broom)
library(tibble)
library(dplyr)
library(tidyr)
library(purrr)

SSlogDistrictCumulativeAdoption100Results$Slopes=1/SSlogDistrictCumulativeAdoption100Results$scale

SSlogDistrictCumulativeAdoption100Results$Intercept=-(SSlogDistrictCumulativeAdoption100Results$xmid/SSlogDistrictCumulativeAdoption100Results$scale)

SSlogDistrictCumulativeAdoption100Results$Year10percent=((-2.2-SSlogDistrictCumulativeAdoption100Results$Intercept)/SSlogDistrictCumulativeAdoption100Results$Slopes)
SSlogDistrictCumulativeAdoption100Results$Year90percent=((2.2-SSlogDistrictCumulativeAdoption100Results$Intercept)/SSlogDistrictCumulativeAdoption100Results$Slopes)
SSlogDistrictCumulativeAdoption100Results$Years10_90=SSlogDistrictCumulativeAdoption100Results$Year90percent-SSlogDistrictCumulativeAdoption100Results$Year10percent

write.csv(SSlogDistrictCumulativeAdoption100Results, file = "SSlogDistrictCumulativeAdoption100Results.csv")

herbi_1st_use <- with(DistrictCumulativeAdoption, seq(2000,2030, length.out = 31))
A.q103_district=levels(as.factor(DistrictCumulativeAdoption$A.q103_district))
herbi_1st_use_dist=expand.grid(herbi_1st_use=herbi_1st_use,A.q103_district=A.q103_district)
herbi_1st_use_dist$pred=predict(SSlogDistrictCumulativeAdoption100, herbi_1st_use_dist)
herbi_1st_use_dist_Data=merge(DistrictCumulativeAdoption,herbi_1st_use_dist,by=c("A.q103_district","herbi_1st_use"),all.y=TRUE)

write.csv(herbi_1st_use_dist_Data, file = "herbi_1st_use_dist_DataPredictionResults100.csv")


Herb_adoption=subset(herbi_1st_use_dist_Data, select=c("A.q103_district","herbi_1st_use","pred") )

names(Herb_adoption)[1:3]=c("Districts","year","Herb_adopt_perc")

write.csv(Herb_adoption, file = "Herb_adoption.csv")



# Plotting predicted herbicide adoption rates

ggplot(Herb_adoption,aes(y=Herb_adopt_perc,x=year,color=Districts)) + 
  geom_line()+
   theme_bw()+
  labs(y="Ädoption rate (% of farmers)")

Code
# Subset the years 2019 to 2030 and three columns 


Herb_adoption_2019_30=subset(Herb_adoption,Herb_adoption$year%in%c(2019:2030))

ggplot(Herb_adoption_2019_30,aes(y=Herb_adopt_perc,x=year,color=Districts)) + 
  geom_line()

Code
Herb_adoption_2019_30$Districts=toupper(Herb_adoption_2019_30$Districts)

APY_plus_adoption=merge(Herb_adoption_2019_30,APY_data_wheat_2019_30,by=c("Districts","year"))


APY_plus_adoption$Area_Hectare_Herb=(APY_plus_adoption$Herb_adopt_perc/100)*APY_plus_adoption$Area_Hectare

2.3 Step 3: Supply, demand elasticities, and prices

Code
APY_plus_adoption$Supply_Elast=0.22 # Taken from Kumar et al 2016
APY_plus_adoption$Demand_Elast=0.340 # Taken from Kumar et al 2011

Prices=import("tabpricedist.csv")

Prices=subset(Prices,select=c("A.q103_district","mean"))

names(Prices)[1:2]=c("Districts","Prices")

Prices$Districts=toupper(Prices$Districts)

APY_plus_adoption=merge(APY_plus_adoption,Prices,by="Districts")

2.4 Step 4: Yield and cost changes due to research

Code
library(rio)
Yield_changes=import("tau.hat_weeding_herb_predictions_dist.csv")

Yield_changes=subset(Yield_changes,select=c("A.q103_district","mean"))

names(Yield_changes)[1:2]=c("Districts","Yield_gain_t_ha")

Yield_changes$Districts=toupper(Yield_changes$Districts)

APY_plus_adoption=merge(APY_plus_adoption,Yield_changes,by="Districts")

# cost changes
APY_plus_adoption$Cost_changes=0

2.5 Step 5: Research induced supply shift

Code
# K_shift_yld_compt

APY_plus_adoption$K_shift_yld_compt=APY_plus_adoption$Yield_gain_t_ha/APY_plus_adoption$Supply_Elast


# K_shift_cost_component
APY_plus_adoption$K_shift_cost_compt=APY_plus_adoption$Cost_changes/(1+(APY_plus_adoption$Yield_gain_t_ha/APY_plus_adoption$Yield_Tonnes_Hectare))

# K_shift
APY_plus_adoption$K_shift=APY_plus_adoption$K_shift_yld_compt-APY_plus_adoption$K_shift_cost_compt

APY_plus_adoption$K_shift_A=APY_plus_adoption$K_shift*APY_plus_adoption$Herb_adopt_perc*0.01

# Z_t=K_t e/e+eta

APY_plus_adoption$Z_t=(APY_plus_adoption$K_shift_A*APY_plus_adoption$Supply_Elast)/(APY_plus_adoption$Supply_Elast+APY_plus_adoption$Demand_Elast)

2.6 Step 6: Compute producer surplus, consumer surplus and total surplus

For a small closed economy, the equations producer surplus, consumer surplus and total surplus.

Producer surplus

\[ \Delta PS_t = P_0Q_0 (K-Z)(1+0.5Z \eta) \]

Consumer surplus \[ \Delta CS_t = P_0Q_0 Z(1+0.5Z \eta) \]

Total surplus \[ \Delta TS_t =\Delta PS_t + \Delta CS_t \]

Code
# Producer surplus
#APY_plus_adoption$producer_surplus=APY_plus_adoption$Production_Tonnes*APY_plus_adoption$Prices*APY_plus_adoption$K_shift_A*(1-0.5*APY_plus_adoption$K_shift_A*APY_plus_adoption$Supply_Elast)

APY_plus_adoption$producer_surplus=APY_plus_adoption$Production_Tonnes*APY_plus_adoption$Prices*(APY_plus_adoption$K_shift_A-APY_plus_adoption$Z_t)*(1+0.5*APY_plus_adoption$Z_t*APY_plus_adoption$Demand_Elast)

APY_plus_adoption$producer_surplus_Rs_per_ha=(APY_plus_adoption$producer_surplus/APY_plus_adoption$Area_Hectare)*1000

# Consumer surplus
APY_plus_adoption$consumer_surplus=APY_plus_adoption$Production_Tonnes*APY_plus_adoption$Prices*APY_plus_adoption$Z_t*(1+0.5*APY_plus_adoption$Z_t*APY_plus_adoption$Demand_Elast)

# Total surplus
APY_plus_adoption$total_surplus=APY_plus_adoption$producer_surplus+APY_plus_adoption$consumer_surplus

2.6.1 Small open economy producer surplus

For a small open economy in which the state’s contribution cannot affect global prices, producer surplus can be computed as: \[ \Delta PS_t = P_0Q_0K_t(1-0.5K_t \epsilon) \]

Code
APY_plus_adoption$producer_surplus_open_econ <- APY_plus_adoption$Production_Tonnes * APY_plus_adoption$Prices * (1 - 0.5 * APY_plus_adoption$K_shift_A * APY_plus_adoption$Supply_Elast)

APY_plus_adoption$producer_surplus_open_econ_Rs_per_ha <- (APY_plus_adoption$producer_surplus_open_econ / APY_plus_adoption$Area_Hectare) * 1000

2.7 Step 7: Cost of research

Code
APY_plus_adoption$Cost_of_research=0

2.8 Step 8: Compute economic indicators

Code
# Discount rate
APY_plus_adoption$discount_rate <- 0.2

# time
APY_plus_adoption$time <- APY_plus_adoption$year - 2019

# cashflows
APY_plus_adoption$net_producer_surplus <- APY_plus_adoption$producer_surplus - APY_plus_adoption$Cost_of_research

APY_plus_adoption$net_producer_surplus_open_econ <- APY_plus_adoption$producer_surplus_open_econ - APY_plus_adoption$Cost_of_research

APY_plus_adoption$net_consumer_surplus <- APY_plus_adoption$consumer_surplus - APY_plus_adoption$Cost_of_research

APY_plus_adoption$net_total_surplus <- APY_plus_adoption$total_surplus - APY_plus_adoption$Cost_of_research

# Calculate the Present Value of Each Cash Flow
APY_plus_adoption$Producer_present_value <- APY_plus_adoption$net_producer_surplus / ((1 + APY_plus_adoption$discount_rate)^APY_plus_adoption$time)

APY_plus_adoption$Producer_present_value_open_econ <- APY_plus_adoption$net_producer_surplus_open_econ / ((1 + APY_plus_adoption$discount_rate)^APY_plus_adoption$time)

APY_plus_adoption$Consumer_present_value <- APY_plus_adoption$net_consumer_surplus / ((1 + APY_plus_adoption$discount_rate)^APY_plus_adoption$time)

APY_plus_adoption$Total_present_value <- APY_plus_adoption$net_total_surplus / ((1 + APY_plus_adoption$discount_rate)^APY_plus_adoption$time)


# Sum all the Present Values
library(data.table)
APY_plus_adoption <- data.table(APY_plus_adoption)

APY_plus_adoption_indic <- APY_plus_adoption[, .(
    Producer_net_present_value_sum = base::sum(Producer_present_value, na.rm = TRUE),
    Producer_present_value_open_econ_sum = base::sum(Producer_present_value_open_econ, na.rm = TRUE),
    Consumer_net_present_value_sum = base::sum(Consumer_present_value, na.rm = TRUE),
    Total_net_present_value_sum = base::sum(Total_present_value, na.rm = TRUE),
    Wheat_area_mean = base::mean(Area_Hectare, na.rm = TRUE),
    Area_Hectare_Herb_mean = base::mean(Area_Hectare_Herb, na.rm = TRUE),
    Production_Tonnes_mean = base::mean(Production_Tonnes, na.rm = TRUE),
    Yield_Tonnes_Hectare_mean = base::mean(Yield_Tonnes_Hectare, na.rm = TRUE)
), by = .(Districts)]

# Producer surplus small open economy
APY_plus_adoption_indic$Producer_present_value_open_econ_sum_annual <- APY_plus_adoption_indic$Producer_present_value_open_econ_sum / 11
sum(APY_plus_adoption_indic$Producer_present_value_open_econ_sum_annual)
[1] 23862234
Code
summary(APY_plus_adoption_indic$Producer_present_value_open_econ_sum_annual)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 201921  452567  776113  917778 1294647 2691572 
Code
APY_plus_adoption_indic$Producer_present_value_open_econ_sum_annual_Rs_ha <- (APY_plus_adoption_indic$Producer_present_value_open_econ_sum_annual / APY_plus_adoption_indic$Wheat_area_mean) * 1000

summary(APY_plus_adoption_indic$Producer_present_value_open_econ_sum_annual_Rs_ha)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  11815   14875   16345   16521   18574   20293 
Code
# Producer surplus
APY_plus_adoption_indic$Producer_net_present_value_sum_annual <- APY_plus_adoption_indic$Producer_net_present_value_sum / 11

sum(APY_plus_adoption_indic$Producer_net_present_value_sum_annual)
[1] 23572896
Code
summary(APY_plus_adoption_indic$Producer_net_present_value_sum_annual)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 196220  451067  784058  906650 1280206 2393513 
Code
APY_plus_adoption_indic$Producer_net_present_value_sum_annual_Rs_ha <- (APY_plus_adoption_indic$Producer_net_present_value_sum_annual / APY_plus_adoption_indic$Wheat_area_mean) * 1000

summary(APY_plus_adoption_indic$Producer_net_present_value_sum_annual_Rs_ha)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  10250   14042   16339   16911   18348   28955 
Code
# Consumer surplus
APY_plus_adoption_indic$Consumer_net_present_value_sum_annual <- APY_plus_adoption_indic$Consumer_net_present_value_sum / 11

sum(APY_plus_adoption_indic$Consumer_net_present_value_sum_annual)
[1] 15253051
Code
summary(APY_plus_adoption_indic$Consumer_net_present_value_sum_annual)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 126966  291867  507332  586656  828369 1548744 
Code
APY_plus_adoption_indic$Consumer_net_present_value_sum_annual_Rs_ha <- (APY_plus_adoption_indic$Consumer_net_present_value_sum_annual / APY_plus_adoption_indic$Wheat_area_mean) * 1000

summary(APY_plus_adoption_indic$Consumer_net_present_value_sum_annual_Rs_ha)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   6632    9086   10572   10943   11872   18735 
Code
# Total surplus
APY_plus_adoption_indic$Total_net_present_value_sum_annual <- APY_plus_adoption_indic$Total_net_present_value_sum / 11

sum(APY_plus_adoption_indic$Total_net_present_value_sum_annual)
[1] 38825947
Code
summary(APY_plus_adoption_indic$Total_net_present_value_sum_annual)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 323186  742934 1291390 1493306 2108575 3942257 
Code
APY_plus_adoption_indic$Total_net_present_value_sum_annual_Rs_ha <- (APY_plus_adoption_indic$Total_net_present_value_sum_annual / APY_plus_adoption_indic$Wheat_area_mean) * 1000

summary(APY_plus_adoption_indic$Total_net_present_value_sum_annual_Rs_ha)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  16882   23128   26911   27854   30220   47690 
Code
# Benefit-cost ratio
# APY_plus_adoption$BCR=APY_plus_adoption$producer_surplus/APY_plus_adoption$Cost_of_research

# Internal Rate of Return (IRR)

# IRR(cf0,cf,times,plot=FALSE)
# MIRR=

# irr1 <- project1_cf %>%
#   select(cf) %>%
#   .[[1]] %>%
#   irr()

# Modified Internal Rate of Return (MIRR)

3 Spatial visualization of the economic surpluses

Code
library(geodata)

India <- gadm(country = "IND", level = 2, path = "shp")
plot(India)

Code
India_aoi <- subset(India, India$NAME_1 == "Bihar")

plot(India_aoi)

plot(India_aoi, add = TRUE)

Code
library(sf)

India_aoi_sf <- st_as_sf(India_aoi)
library(mapview)

mapview(India_aoi_sf)
Code
India_aoi_sf$Districts <- toupper(India_aoi_sf$NAME_2)

APY_plus_adoption_indic_sf <- merge(India_aoi_sf, APY_plus_adoption_indic, by = "Districts")

# Average wheat area
# Using tmap
library(tmap)
tmap_mode("plot")

Wheat_area_sf_tmap <- tm_shape(APY_plus_adoption_indic_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Wheat_area_mean", palette = "viridis", scale = 2.5, title = "Wheat area (ha)")
Wheat_area_sf_tmap

Code
tmap_save(Wheat_area_sf_tmap, "figures/Wheat_area_sf_tmap.png", dpi = 300, width = 4.88, height = 3.16)

# Average wheat yields
tmap_mode("plot")

Wheat_yield_sf_tmap <- tm_shape(APY_plus_adoption_indic_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Yield_Tonnes_Hectare_mean", palette = "viridis", scale = 2.5, title = "Wheat yields (t/ha)")
Wheat_yield_sf_tmap

Code
tmap_save(Wheat_yield_sf_tmap, "figures/Wheat_yield_sf_tmap.png", dpi = 300, width = 4.88, height = 3.16)


# Average wheat production
tmap_mode("plot")

Wheat_production_sf_tmap <- tm_shape(APY_plus_adoption_indic_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Production_Tonnes_mean", palette = "viridis", scale = 2.5, title = "Wheat production (t)")
Wheat_production_sf_tmap

Code
tmap_save(Wheat_production_sf_tmap, "figures/Wheat_production_sf_tmap.png", dpi = 300, width = 4.88, height = 3.16)



# Producer surplus

mapview(APY_plus_adoption_indic_sf, zcol = "Producer_net_present_value_sum_annual", layer.name = "Producer surplus (INR 1000)")
Code
# Using tmap
library(tmap)
tmap_mode("plot")

prod_surplus_changes_sf_tmap <- tm_shape(APY_plus_adoption_indic_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Producer_net_present_value_sum_annual", palette = "viridis", scale = 2.5, title = "Producer surplus (INR 1000)")
prod_surplus_changes_sf_tmap

Code
tmap_save(prod_surplus_changes_sf_tmap, "figures/prod_surplus_changes_sf_tmap.png", dpi = 300, width = 3.88, height = 3.16)

# Small open economy
library(tmap)
tmap_mode("plot")

prod_surplus_changes_open_econ_sf_tmap <- tm_shape(APY_plus_adoption_indic_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Producer_present_value_open_econ_sum_annual", palette = "viridis", scale = 2.5, title = "Small open economy PS (INR 1000)")
prod_surplus_changes_open_econ_sf_tmap

Code
tmap_save(prod_surplus_changes_open_econ_sf_tmap, "figures/prod_surplus_changes_open_econ_sf_tmap.png", dpi = 300, width = 3.88, height = 3.16)



# Producer surplus per ha
# Using tmap
library(tmap)
tmap_mode("plot")

prod_surplus_changes_perha_sf_tmap <- tm_shape(APY_plus_adoption_indic_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Producer_net_present_value_sum_annual_Rs_ha", palette = "viridis", scale = 2.5, title = "Producer surplus per ha")
prod_surplus_changes_perha_sf_tmap

Code
tmap_save(prod_surplus_changes_perha_sf_tmap, "figures/prod_surplus_changes_perha_sf_tmap.png", dpi = 300, width = 3.88, height = 3.16)

## Small open economy
library(tmap)
tmap_mode("plot")

prod_surplus_changes_open_econ_perha_sf_tmap <- tm_shape(APY_plus_adoption_indic_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Producer_present_value_open_econ_sum_annual_Rs_ha", palette = "viridis", scale = 2.5, title = "Small open economy PS per ha")
prod_surplus_changes_open_econ_perha_sf_tmap

Code
tmap_save(prod_surplus_changes_open_econ_perha_sf_tmap, "figures/prod_surplus_changes_open_econ_perha_sf_tmap.png", dpi = 300, width = 3.88, height = 3.16)



mapview(APY_plus_adoption_indic_sf, zcol = "Consumer_net_present_value_sum_annual", layer.name = "Consumer surplus (INR 1000)")
Code
mapview(APY_plus_adoption_indic_sf, zcol = "Total_net_present_value_sum_annual", layer.name = "Total economic surplus (INR 1000)")
Code
library(reactable)
library(htmltools)
library(fontawesome)

htmltools::browsable(
    tagList(
        tags$button(
            tagList(fontawesome::fa("download"), "Download as CSV"),
            onclick = "Reactable.downloadDataCSV('APY_plus_adoption_indic_sf', 'APY_plus_adoption_indic_sf.csv')"
        ),
        reactable(
            APY_plus_adoption_indic_sf,
            searchable = TRUE,
            defaultPageSize = 38,
            elementId = "APY_plus_adoption_indic_sf"
        )
    )
)
Code
# Yield changes
Yield_changes <- subset(Yield_changes, !(Yield_changes$Districts == "KISHANGANJ"))

library(mapview)
Yield_changes_sf <- merge(India_aoi_sf, Yield_changes, by = "Districts")

mapview(Yield_changes_sf, zcol = "Yield_gain_t_ha", layer.name = "Yield gain to herbicides (tons per ha)")
Code
# Using tmap
library(tmap)
tmap_mode("plot")

Yield_changes_sf_tmap <- tm_shape(Yield_changes_sf) +
    tm_borders(alpha = 0.5, col = "white") +
    tm_polygons(col = "Yield_gain_t_ha", palette = "viridis", scale = 2.5, title = "Yield gain to herbicides (tons per ha)")
Yield_changes_sf_tmap

Code
tmap_save(Yield_changes_sf_tmap, "figures/Yield_changes_sf_tmap.png", dpi = 300, width = 3.88, height = 3.16)



library(reactable)
library(htmltools)
library(fontawesome)

htmltools::browsable(
    tagList(
        tags$button(
            tagList(fontawesome::fa("download"), "Download as CSV"),
            onclick = "Reactable.downloadDataCSV('Yield_changes_sf', 'Yield_changes_sf.csv')"
        ),
        reactable(
            Yield_changes_sf,
            searchable = TRUE,
            defaultPageSize = 38,
            elementId = "Yield_changes_sf"
        )
    )
)

4 Conclusion

The notebook has demonstrated the workflow for calculating spatially disaggregated economic surplus measures that can be used to target and prioritize agricultural interventions.

5 References

Alston, J.M., Norton, G.W., and Pardey, P.G. 1995. “Science under scarcity: Principles and practice for agricultural research evaluation and priority setting.” CAB International. Wallingford, UK.